Import packages

rm(list = ls())
library(maps)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(mapdata)
library(ggplot2)

source("../R/geomatching.R")

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

source("../R/idw2hr.R")

## Loading required package: iterators

## Loading required package: parallel

source("../R/hr2poly.R")

Import data

• AQ_EEA_NO2: NO₂ point data measured from air quality monitoring data.

• AQ_CAMS_N02: NO₂ gridded data (0.01° x 0.01°) obtained from CAMS mathematical model.

# CSV data (.csv)
AQ_EEA_NO2 <- read.csv("AQ_EEA_NO2.csv")

# R environment (.RData)
load("AQ_CAMS_NO2.RData")

1^st step: geomatching

Configuration

data_dict <- list(
  "AQ_EEA_NO2" = list(
    "data" = AQ_EEA_NO2,
    "format" = "xyt",
    "type" = "points",
    "crs" = 4979
  ),
  "AQ_CAMS_NO2" = list(
    "data" = AQ_CAMS_NO2,
    "format" = "matrix",
    "type" = "grid",
    "crs" = 4326
  )
)

data_ls <- list()
settings = list(
  "format"=c(),
  "type"=c(),
  "crs"=c()
)
for (i in seq_along(data_dict)) {
  data_ls[[i]] <- data_dict[[i]]$data
  settings$format[[i]] <- data_dict[[i]]$format
  settings$type[[i]] <- data_dict[[i]]$type
  settings$crs[[i]] <- data_dict[[i]]$crs
}

Run

res_geomatch <- geomatching(
  data=data_ls,
  settings = list(
    "format"=settings$format,
    "type"=settings$type,
    "crs"=settings$crs
  )
)

## [1] "reading data 1"
## [1] "done"
## [1] "reading data 2"
## [1] "done"
## [1] "converting data 1 to ST"
## [1] "done"
## [1] "converting data 2 to ST"
## [1] "done"

# rename columns
res_geomatch <- res_geomatch %>%
  rename(
    altitude = Altitude,
    AQ_CAMS_NO2 = matrix_2
  )

Data visualization

for (var in c("AQ_EEA_NO2", "AQ_CAMS_NO2", "altitude")) {
  
  p <- ggplot() +
    geom_map(
      data = map_data("italy"),
      map = map_data("italy"),
      aes(map_id = region),
      fill = "lightgrey",
      color = "black",
      linewidth = .1,
      alpha=.9
    ) + 
    geom_point(
      data = res_geomatch[res_geomatch$time==res_geomatch$time[1],],
      aes(x = longitude, y = latitude, color = .data[[var]]),
      size = 5,
      
    ) +
    labs(
      title = sprintf("Geographical map of %s", var),
      subtitle = "Output of geomatching",
      x = "Longitude",
      y = "Latitude"
    ) + 
    scale_color_viridis_c(
      option = "C"
    ) +
    theme(
      legend.title = element_blank(),
      plot.title = element_text(face = "bold")
    ) 
  print(p)
}

2^nd step: idw2hr

res_idw2hr <- idw2hr(
  data = res_geomatch,
  col_names = list(
    "lon" = "longitude",
    "lat" = "latitude",
    "t" = "time"
  ),
  interest_vars = list("altitude", "AQ_EEA_NO2", "AQ_CAMS_NO2"),
  crs = 4979
)

Data visualization

for (var in c("AQ_EEA_NO2", "AQ_CAMS_NO2", "altitude")) {
  
  p <- ggplot() +
    geom_point(
      data = res_idw2hr[res_idw2hr$time==res_idw2hr$time[1],],
      aes(x = longitude, y = latitude, color = .data[[var]]),
      size = 5,
      
    ) +
    geom_map(
      data = map_data("italy"),
      map = map_data("italy"),
      aes(map_id = region),
      fill = "lightgrey",
      color = "black",
      linewidth = .1,
      alpha=.05
    ) +
    labs(
      title = sprintf("Geographical map of %s", var),
      subtitle = "Output of idw2hr",
      x = "Longitude",
      y = "Latitude"
    ) + 
    scale_color_viridis_c(
      option = "C"
    ) +
    theme(
      legend.title = element_blank(),
      plot.title = element_text(face = "bold")
    ) 
  print(p)
}

3^rd step: hr2poly

res_hr2poly <- hr2poly(
  res_idw2hr,
  polygon_type = "mun",
  crs = 4979
)

Data visualization

for (var in c("AQ_EEA_NO2_min", "AQ_EEA_NO2_mean", "AQ_EEA_NO2_median",
              "AQ_EEA_NO2_max", "AQ_EEA_NO2_std", "AQ_CAMS_NO2_min", "AQ_CAMS_NO2_mean",
              "AQ_CAMS_NO2_median", "AQ_CAMS_NO2_max", "AQ_CAMS_NO2_std", "altitude_min", 
              "altitude_mean", "altitude_median", "altitude_max", "altitude_std")) {
  
  p <- ggplot() +
    geom_sf(
      data=res_hr2poly,
      aes_string(fill=var)
    ) +
    labs(
      title = sprintf("Geographical map of %s", var),
      subtitle = "Output of hr2poly",
      x = "Longitude",
      y = "Latitude"
    ) + 
    scale_fill_viridis_c(
      option = "C"
    ) +
    theme(
      legend.title = element_blank(),
      plot.title = element_text(face = "bold")
    )
  print(p)
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.